%% Analyze for thresholds based on reversals
Nrevs = 10; % number of reversals

rp = schedule.readparams;
rv = schedule.response_vals;

if isnan(rv(end)), rv(end,:) = []; end

ttypidx = findincell(strfind(rp,'trial_type'))+1;
codeidx = findincell(strfind(rp,'response_code'))+1;
levlidx = findincell(strfind(rp,'ToneDB'))+1;
freqidx = findincell(strfind(rp,'ToneFreq'))+1;

ufreq = unique(rv(:,freqidx));

fprintf('\nFreq (Hz)\tThresh (dB)\tFA (%%)\t#Reversals\n')
clear data
for i = 1:length(ufreq)   
    ind = rv(:,freqidx) == ufreq(i);
    
    clear d
    
    % trial levels
    d.trials.level = rv(ind,levlidx);
    d.trials.ind   = ind;
    d.trials.catch = ind & rv(:,ttypidx);
    ctind = ind&rv(:,ttypidx);
    d.FArate = sum(bitget(rv(ctind,codeidx),3))/sum(ctind)*100;
%     d.trials.lwocind  = ind & ~rv(:,ttypidx);
%     d.trials.lwoc  = rv(d.trials.lwocind,levlidx);
    
    % reversals at this frequency
    k = 1;
    level = d.trials.level;
    if numel(level) == 1
        d.reversal.level = level;
        d.reversal.idx = 1;
        d.reversal.direction = nan;
    else
        d.reversal.level = nan;
        d.reversal.idx = nan;
        d.reversal.direction = nan;
        lastdir = level(2) > level(1); % down: 0 , up: 1
        for j = 2:length(level)
            if lastdir && level(j) < level(j-1) ...
                    || ~lastdir && level(j) > level(j-1)
                d.reversal.level(k) = level(j-1);
                d.reversal.idx(k) = j-1;
                d.reversal.direction(k) = lastdir;
                lastdir = ~lastdir;
                k = k + 1;
            end
        end
        % last level is always a reversal
        d.reversal.level(end+1) = level(end);
        d.reversal.idx(end+1)   = length(level);        
    end
    
    if length(d.reversal.level) < Nrevs
        d.threshold = nan;
        d.reversal.residual = nan;
    else
        d.threshold = mean(d.reversal.level(end-Nrevs+1:end));
        d.reversal.residual = d.reversal.level - d.threshold;
    end
        
    d.code = rv(ind,codeidx);
    d.freq = ufreq(i);
    fprintf('% 0.0f\t\t\t%0.0f\t\t\t%0.0f\t\t\t%d\n', ...
        ufreq(i),d.threshold,d.FArate,numel(d.reversal.level))
    data(i,1) = d; %#ok<SAGROW>
end
thr = [data.threshold];
clipboard('copy',sprintf('%d\n',thr))
fprintf('Data copied to clipboard\n\n')


%% Plot
clf
for i = 1:length(data)
    subplot(6,2,i);
    plot(data(i).trials.level,'-ob','markerfacecolor','b')
    axis tight
    ylim([-20 70]);
    
    hold on
    plot(data(i).reversal.idx,data(i).reversal.level,'ob', ...
        'markerfacecolor','g');
    plot(xlim,[data(i).threshold data(i).threshold],':k');
    hold off
    
    title(sprintf('%0.1f kHz, %0.0f dB (%0.0f)', ...
        data(i).freq/1000,data(i).threshold,data(i).FArate));
end
subplot(6,2,9:12)
plot([data.freq]/1000,[data.threshold],'-ok','markerfacecolor','k','markersize',10, ...
    'linewidth',3);
xlabel('Frequency (kHz)'); ylabel('Level (dB SPL)');
title('Audiogram')
ylim([-20 70]);
hold on
plot([data.freq]/1000,[data.FArate],'*r');
plot(repmat([data.freq]/1000,2,1),[zeros(1,length(data)); [data.threshold]],'-k');
plot(xlim,[0 0],'-k');
hold off

